Apparent Fracture in Polymeric Fluids under Step Shear 
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Recent step strain experiments in well-entangled polymeric liquids demonstrated a bulk fracture- 
like phenomenon. This was taken, in these experiments, to signify the invalidity of the Doi-Edwards 
(DE) tube model. We have investigated this phenomenon using the Rolie-Poly equation, which 
approximates a successful version of the DE theory, and we find close quantitative agreement with 
the experiments, as well as with the proposal by Marrucci and Grizzuti in 1983 that entangled 
polymer liquids possess an elastic instability. The fracture is a transient manifestation of this elastic 
instability that relies on the amplification of spatially inhomogeneous fluctuations. It resembles 
spinodal decomposition, with strain playing the role of the conserved quantity. 
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Introduction — Entangled polymers are ubiquitious 
in the plastics industry, and control of their processing is 
vital for producing well-tailored materials [1]. The mo- 
tion of an entangled polymer molecule is restricted to a 
tube-like region due to the constraints imposed by sur- 
rounding chains. The theory for this, due to Doi and Ed- 
wards (DE)[2], predicts a maximum in the shear stress 
T xy as a function of shear rate [Fig. 1(a)], at a shear rate 
7 roughly equal to the reciprocal of the time Td for a 
polymer to diffuse along its tube (reptation time). This 
non-monotonic constitutive behaviour indicates instabil- 
ity, which can lead to inhomogeneous flows and shear 
banding [3] . This constitutive instability was widely im- 
plicated [4] in the spurt effect [5] , responsible for instabil- 
ities in industrial processes; however, spurt is now usu- 
ally attributed to wall slip [6]. In rapid startup flow the 
DE theory also predicts the rubber elastic behaviour of 
a stress overshoot [2, 7]. 

Banding was not inferred in early shearing experiments 
on polymer melts [8], and subsequently the mechanisms 
of chain stretch and convected constraint release (CCR) 
- chain relaxation due to the release of entanglement con- 
straints - were incorporated into the DE theory [9]; CCR 
can restore stable constitutive behavior. However, new 
observations of shear banding in polymer solutions [10] 
seem to validate the DE instability [7, 11], and it remains 
unknown how active CCR is. 

Our motivation is a surprising observation by Boukany 
et al. in step strain experiments (using flat plates) in 
poly(styrene-butadiene) melts with Z ps 53 — 160 chain 
entanglements and strains 70 > 2 [12], at very fast shear 
rates. After the step they observed homogeneous relax- 
ation for a short time, after which the stress relaxed more 
rapidly and the material split into two layers moving in 
opposite directions, separated by a thin [< 40 /jm] 'frac- 
ture' layer [Fig. 1 of [12]]. They claim that the tube 
model cannot describe this behaviour [7, 11]. A similar 
observation was reported in poly(ethylene oxide) melts 
[13]. However, long ago Marrucci and Grizzuti (MG) 
showed that the DE theory predicts non-uniform elastic 



response for large enough rapid step strains [7], which 
could explain the anomalously fast stress relaxation of 
some melts subjected to step strains [14, 15]. Here we 
show how fracture follows from the elastic nature of the 
DE instability [7]. 

Model — We separate the total stress tensor T into 
contributions from the polymer and a Newtonian solvent, 



GW + r/(n + k t ) -pi, 
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where G is the plateau modulus, 77 is the solvent viscos- 
ity, the pressure p maintains incompressibility, I is the 
identity tensor and K a p = dv a /drp. The fluid velocity v 
obeys 
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where p is the fluid density. For very small Reynolds 
number, Eq. (2) reduces to V ■ T = 0, which is used for 
the spatially resolved calculations. 

The dimensionless polymeric conformation, or strain, 
tensor W is assumed to obey the diffusive Rolie-Poly 
(RP) equation [7, 16], 
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which is a simplified form of the GLaMM model for linear 
polymers [17]. Here, T> is the stress diffusion constant, 
Td is the reptation time, and the Rouse time tr governs 
the relaxation of stretch Tr(W). The parameter /3 quan- 
tifies convective constraint release (CCR); a large value 
of f3 corresponds to more CCR, which leads to mono- 
tonic (stable) behaviour of the shear stress. Following 
successful fits to experiments [16] we use 6 = —0.5. 

Calculations — We consider two infinite flat plates of 
separation L where the top plate is free to move and the 
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FIG. 1. (a) Constitutive (solid line) and steady state shear banding (black dashes) curves. The stress overshoot is indicated 
by green squares (□); the homogenous base state becomes unstable (cj max > 0) during startup at stresses given by the blue 
diamonds (()). (b) Velocity profile at t$ = 0.01249^ (just before shear cessation) for (7) = 200. (c) Stress relaxation for step 
strains 70 = 0.2,2.5; the solid line is for 70 = 2.5 with no initial perturbation. The dot-dash line shows the evolution of the 
most unstable eigenvalue u max , which becomes unstable (o; m ax > 0) in the red (dashed) region, (d) Velocity profiles during 
fracture, with experimental data from [12] superposed, (e) Shear rate profiles, (f) stress relaxation, and (g) evolution of the 
maximum stretch in the gap Tr A max . Parameters: Z = 72, tr = r d /216, (7) = 200, 70 = 2.5, and to = 0.01250r d . Time t 
displayed in units of r^.] 



bottom plate is fixed. The flow direction is x and the 
gradient direction is y so that v = v x (t,y)x and W = 
W(t, y). The following dimensionless quantities are used: 
7 = jr d , V = Vr d jL\ e = r}/(Gr d ), p = {pL 2 )/{Grj), 
v = { T dV) I ' L, and t = t/r d . The degree of entanglement 
Z determines the Rouse time via tr = r d /(3Z) [16-18]. 
A desired average shear rate is imposed for a duration to 
leading to a strain 70 = (7) to- 

The values r d = 310 s and Z = 55 — 100 are consis- 
tent with the data in [12]; with 77 w lPas and G ps 
7 x 10 3 Pa [19] we find e s» 10 -7 ; for numerical stability 
we use e = 10 -4 . For L = 1 mm, p 10 3 kg m~ 3 gives 
p w 10~ 10 and we use V — 10 -5 [20]. Spatial derivatives 
are discretized using a semi implicit central finite differ- 
ence scheme. For a timestep St = 10~ 6 and 1000 spatial 
mesh points the maximum velocity in the fracture con- 
verges within 1 — 4% compared to 5t = 10~ 5 or a mesh 
of 500. The time to fracture convergences similarly. 

We analyze stability using u(t) = [7, A IX , A xy , A yy ] — 
u(t) + J2 k 5u k (t) exp(ifcy), where A = W - I. The ho- 
mogeneous base state u(t) is obtained by solving Eq. (3) 
for T> — 0. Substituting u(t) into Eqs. (2, 3) and lin- 
earizing yields S\ik — M k -Su k . We define uj max as the 
largest real part of the spectrum of eigenvalues of M&; 
for uj max > the base state u(i) is unstable with respect 
to small amplitude fluctuations [21, 22]. 

To capture the behaviour reported in [12], we con- 



sider a fluid with non-monotonic constitutive behaviour, 
f3 = [solid line in Fig. 1(a)], and use Z = 72 (consis- 
tent with [12]); this leads to shear banding and a stress 
plateau in steady state [dashes in Fig. 1(a)] [3]. We ini- 
tialize Eq. (3) with random perturbations of the form 
Su(0,y) = ^Yfn=i( A n/ n2 ) c osmry, A n e [-1,1], where 
£ sets the overall scale of the perturbation. The penalty 
1/n 2 arises because high wavenumbers n should be sup- 
pressed by both spatial gradients in W and by the slow 
dynamics of long wavelength velocity fluctuations that 
induce perturbations upon sample loading (for example) . 
We use £ = 0.01, consistent with the scale of typical ther- 
mal fluctuations in W [23]. 

Perturbations can grow if the fluid becomes unsta- 
ble [7, 22, 23]. Different random perturbations applied 
simultaneously to 7 and A result in different velocity 
profiles after shear cessation at to [24]. For 34% of 300 
such simulations, the resulting velocity profiles were sim- 
ilar to the case reported in [12]. In all cases, the velocity 
profile was nearly homogeneous before shear cessation 
[Fig. 1(b)]. Using initial conditions (see later) that pro- 
duce the experimentally observed velocity profile, we sim- 
ulate examples reported in [12]: (I) Intermediate shear 
rates (7)771 « 1 and an applied strain 70 exceeding the 
strain j ov at which a stress overshoot occurs; (II) high 
shear rates (7)17? > 1 in the stretching regime and moder- 
ate strains 70 < "fov'i and (III) lower shear rates (7)772 < 1 
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FIG. 2. Spatial profiles of (a) initial perturbation; (b) lo- 
cal strain; (c) effective modulus A cfi as well as the unsta- 
ble growth rate u; ma x, after cessation of flow and subsequent 
stretch relaxation; (d) Evolution of unrelaxed polymer seg- 
ments )J,(y,t) during fracture development. [Parameters as in 
Fig. 1. Time t displayed in units of r<j.] 



and 7o < Jov In all cases (7)1^ > 1 and, as shown in 
Fig. 1(a), a homogeneous state is unstable before the ces- 
sation of flow at the imposed strain 70. 

Case I: (7)^ ~ 1 and 70 > 7 „ - We impose (7) = 
200 ((7)17? = 0.93) for a strain 70 = 2.5. Immediately 
before cessation, at t$ , the velocity profile is impercepti- 
bly inhomogeneous [Fig. 1(b)], while at f$ the fluid stops 
with a slight inhomogeneity induced by the perturbation 
[Fig. 1(d)]. This perturbation slowly grows and localizes, 
leading to a 'fracture' plane at which the fluid shears 
very rapidly [Figs. l(de)] and a sizeable stretch Tr A is 
induced [Fig. 1(g)] [24]. 

Some stress quickly relaxes due to stretch relaxation 
in a time t s ~ 7tr [Fig. l(cfg)] just after shear cessation; 
followed by an induction time U ~ 30tr that resembles 
homogenous reptation [blue circles in Figs. 1(c), 1(f)], 
before relaxing quickly in a time t/ — 15tr during the 
'fracture' [Fig. 1(d)]. Thereafter it relaxes like a quies- 
cent melt with a small initial strain 70 = 0.2 [Fig. 1(c)], 



so that the stress released in the fracture corresponds 
to a plastic strain j p \ ~ 2.3. Since the boundaries are 
fixed, positive shear strain within the slip layer is bal- 
anced by opposing recoil in the still-entangled outer re- 
gions [e.g. Fig. 1(e) for t/ra > 0.15]. Without an initial 
perturbation only quiescent relaxation obtains [solid line 
of Fig. 1(c)]. The velocity profiles [Fig. 1(d)] are con- 
sistent with Fig. 1 of [12] (which has an induction time 
k « 5tr). 

Stability — Fig. l(ac) shows that the material is un- 
stable (w max > 0) from well before the stress overshoot 
until shear cessation. To understand this instability, we 
turn to the Marrucci-Grizzuti (MG) observation that, 
for strain 70 > 2.1 the elastic energy function ^(7) for 
the DE model has a negative effective shear modulus 
A = d 2 F/dj 2 < [7], which heralds instability. Thus 
MG predicted elastic instability for a step strain, for 



A ett = u(t + t s ) 



d 2 F 
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where u(t) is the fraction of unrelaxed tube segments, 
whose value decreases when the chain segments lose their 
anisotropy and thus relax stress [2] 

The anisotropy of the polymer conformation tensor W 
defines fj, = \X% — A2 j / |Ai + A2I, where Aj are the eigen- 
values of W in the shear plane [25] . For a homogeneous 
state subjected to a step strain fi(t) relaxes homoge- 
neously to zero. However, an inhomogeneous initial con- 
dition initiates instability and an inhomogeneous fi(y, t) 
[Fig. 2(d)], such that the fracture region experiences en- 
hanced anisotropy and polymer stretch [Fig. 1(g)]. 

Fig. 2 (be) shows the spatial profiles for the strain and 
the effective shear modulus A after stretch relaxation 
[26]. The fracture region is most unstable, so that the 
initial perturbation [Fig. 2(a)] can localize strain. In the 
elastic limit rr d > 1, A cS ~ dT xy /dj = j^dT^/dt < 
[7, 11, 22, 23], which coincides with the stress overshoot. 
The unstable region predicted by the elastic limit coin- 
cides with the most unstable eigenvalue w max calculated 
from the full dynamics, which indicates instability before 
the stress overshoot is reached [e.g. Fig. 1(a)] because of 
the viscous contribution to the instability [22]. 

Conditions for fracture — A detailed study shows 
that perturbations in A xx and A yy induce fracture [24]. 
The step strain 70 advects the initial perturbation into 
into a shear component of the polymer strain [e.g. 
W xy (y,t ) ^ 7o(l + Ayy{y,ty), which generates an in- 
homogeneous shear rate #7(2/,^) ~ —j o A yy (y,0)/e im- 
mediately after cessation of flow to maintain V • T ~ 0]. 
Although general perturbations are complex [Fig. 2(a)] 
[24] , if the perturbation leads to a local maximum in the 
polymeric strain 7, this then defines the position with 
most negative effective shear modulus A e g < [Fig. 2(c)] 
[26], and thus the fracture position. This is born out by 
the dynamical calculation [Fig. 2(c)], which shows a criti- 
cal unstable eigenvector dominated by the growth of A xx 
and hence enhanced stretch in the flow direction [24]. 
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FIG. 3. (a) Values of Z and 70 required for fracture at fixed (7); contours show the maximum shear rates during fracture, (be) 
Case II ((7} = 900, (7)7^ = 4.2, 70 = 2.5): (b) Stress relaxation and unstable growth rate U) m&K ; (c) velocity profiles during 
fracture, (d) Stress decay for three different imposed strains 70. (ef) Case III ((7) = 10, (7)7^ = 0.046, 70 = 1.3): (e) Stress 
relaxation and u max , and (f) velocity profiles. [All other parameters as in Fig. 1. Time t displayed in units of r c ;.] 



The subsequent evolution resembles spinodal decom- 
position of a conserved quantity, since the total strain j 
is fixed. The strain in the most unstable region grows 
while that in the less unstable regions decreases, thus 
reducing their degree of instability. This leads to recoil 
and a sharpening of the deformation around the most 
unstable position, which can then fracture. The growing 
local maximum competes with reptation which relaxes 
the instability; hence for fracture to occur, the initial 
amplitude must be sufficiently large for it to localize and 
induce fracture before reptation relaxes the system. Flu- 
ids with large amounts of convected constraint release do 
not undergo fracture, because of the enhanced relaxation. 

Character of Fracture — A larger strain leads to a 
less dramatic fracture [Fig. 3(ad)] because the total stress 
has passed the overshoot and decreased, hence releasing 
less stress into the fracture; however the molecular strain 
W xy is larger, which leads to a faster growing instability. 
This is consistent with Fig. 8 of [12]. Alternatively, for 
a higher imposed strain rate the response is more domi- 
nated by stretching, which leaves less orientational stress 
and molecular strain after stretch relaxation so that frac- 
ture that takes longer to develop [22]. 

Further Comparisons: - We present two more com- 
parisons with Ref. [12]. In Case II ((7) = 900, (7)7* 
= 4.2), the shear rate is large but the strain 70 = 2.5 
is slightly less that the overshoot strain 7 „ [Fig. 3(bc)]. 
The velocity profiles are similar to Case I, and consistent 
with Fig. 2 of [12]. Because the growth rate w max is so 
rapid for the high shear rate, the smaller strain can effect 
the necessary large growth of the instability. In this case 
the induction time t t ~ SQtjj is similar to that of Case I. 

In Case III ((7) = 10, (7)7^ = 0.046) the shear rate 



is relatively weak and the applied strain is just less than 
that of the stress overshoot 70 < jov [Fig. 3(ef)]. Here, 
the 'fracture' and recoil is much weaker than for the pre- 
vious two cases, due to the small growth rate. Moreover, 
the stress does not show any distinguishable delay since 
the shear rate is so low that the fluid response is vis- 
cous dominated and its elastic nature does not show up 
clearly. The weak recoil agrees with Fig. 7 of [12]. 

Fig. 6 of [12] demonstrated that, for sub-overshoot 
strains, higher shear rates lead to longer induction times. 
However, our calculations predict that higher shear rates 
lead to longer induction times because of the faster grow- 
ing instability [24]. We cannot explain this discrepancy. 

Conclusion — Our calculations suggest that the frac- 
ture seen in recent step strain experiments on poly- 
meric liquids [12, 13] results from the underlying elas- 
tic instability whose signature is stress overshoot during 
rapid startup [7, 27]. Once stretch degrees of freedom 
have relaxed, the deformed melt is elastically unstable 
so that small inhomogeneities grow into plastic strain 
(shear flow) in the most unstable regions. If this insta- 
bility grows fast enough compared to reptation then a 
dramatic fracture can result for particular shapes of the 
initial disturbance. The perturbation's shape and ampli- 
tude control whether fracture occurs. 

In related work, Manning et al. studied a similar 
phenomenon in a shear-transformation-zone model of an 
amorphous solid [28], which can yield plastically via a 
fluid 'shear band' during startup of steady shear flow. 

Boukany et al. suggested that the fracture demands 
new physics [12]. Certainly current tube models are in- 
complete [29]. However, our calculations are reasonable 
if the spatial features are smooth on length scales greater 
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than the tube diameter a ~ 3 — 4nm. For a gap of 1mm, scription of the tube model, as used here, is adequate, 
the fracture width 8x ~ 0.05 corresponds to a thickness 

of order 50 fim, which is consistent with the dimension Acknowledgments — This study was funded by the 
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Apparent Fracture in Polymeric Fluids under Step Shear: Supplementary Information 

I. CALCULATIONS 
A. Step-Strain Calculations for Different Initial Conditions 

The starting point for the calculations is the diffusive Rolie-Poly (DRP) model, given by [7, 16] 



^ = ..w + w.^-l(w-i)- 2 ^^^ fw + if^Vcw-Di+pyw, (5) 

at r d t r , y \ 3 J J 

where W is a polymer strain, K a p — d a vp, v is the fluid velocity, and tr are the reptation and stretch relaxation 
times respectively, I is the identity tensor, (3 measures the amount of convective constraint release in the system, 8 is 
a fitting parameter and T> is stress diffusion constant. 

We use the Cartesian coordinate system (for the case of simple shear flow where the fluid is placed between 
two infinite parallel plates of separation L) where y is the velocity gradient direction and x is the flow direction, 
v = v x (t, y)H and W = W(i, y). Substitution into Eq. 5 with W = A + I gives 



<>Xj 2A xy i - A„- — [1 - A] [(PA + l)A xx + l]+V^ A 
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8t xyl xx T R l JL ^ ' ' xx ' J ' dy 

dAxy = $ + ^A yy -A xy -^[1- A]((3A + l)A xy + 5 d2Axy 
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A vy - ^[1 - A] [{fiA + l)A yy + 1}+V- 
TrAV 1/2 



2A (6) 

^yy 

m - -yy TR l^ - JL VM--r^-W-r*j-r- Qy2 



A= 1 



IT) ■ 

where 7 = jrd, T> = T>Td/L 2 . The total stress T is then obtained from W and a Newtonian solvent of viscosity rj as 

T = GW + t}(k + k t ) - pi, (7) 
where G is the plateau modulus and p is pressure, this gives the total shear stress as 

T xy = GA xy + 777. (8) 

To capture the behaviour reported in [12], we initialize Equations 6 with random perturbations of the form 

5 

5u(Q,y) = e^(A„/n 2 )cos(n7ry), (9) 

n=l 

where u = [7, A xx , A xy , A TO ]. The amplitudes A n are chosen randomly within [—1, 1]. The parameter £ = 0.01 sets 
the overall scale of amplitude and a cosine series was chosen since it satisfies the boundary condition imposed on A. 
Using more modes does not change the resultant perturbation significantly due to the 1/n 2 penalty on the amplitudes. 
Each component of u is initially perturbed separately using different random perturbations. Then all components of 
u are perturbed together with each quantity receiving a separate random perturbation. Sample results from these 
simulations are shown in Figs. 4, 5 and C. 
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FIG. 4. Recoil or fracture for different random initial conditions upon perturbing A xx (Left, about 31% fracture) and A yy 
(Right, about 44% fracture). In all cases, the blue line is the perturbation and the red line is the velocity profile when both 
fi v ± reach their extrema together (for fracture) or [i v ± reach their extrema separately (for recoil without fracture). Left axis: 
perturbation; Right axis: velocity. 
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FIG. 5. Perturbation of A xy (Left) and 7 (Right). Only a weak recoil or a weak sign of fracture is seen in both cases. The 
blue and red lines have the same meaning as in Fig. 4. Left axis: perturbation; Right axis: velocity. 



In all calculations reported here, the parameters were set as Z = Td/(3r^) = 72, T> = 1CP 5 , e = rj/(Grd) = 1CP 4 , 
/? = and for stability analysis, p = 10~ 10 . To determine if fracture has occurred or not, consider the 'velocity 
moments' fi v ±, defined by 



i 



(10) 
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FIG. 6. Recoil or fracture upon perturbing all components, with each component receiving a separate random perturbation. 
Red line: perturbation to A„. Green line: perturbation to A yy . Blue line: perturbation to A xy . Magenta line: perturbation 
to 7. Cyan line: recoil or fracture velocity profile v. Left axis: perturbation; Right axis: velocity. The left hand side shows 
about 56% fracture while the right hand side shows about 44% fracture (subfigure P2 shows an example of two fracture planes) . 
In all cases, the velocity profiles occur at the time specified in Fig. 4. 




FIG. 7. Spinodal (shear stress at which the homogeneous state goes unstable), constitutive curve, stress plateau and the 
overshoot stress for the DRP model with /3 = 0, Z = Td/i^Tii) = 72. 



where the sum is over all spatial positions yi and H is the Heaviside step function. If both positive moment /i v+ 
and negative moment fi v _ occur together at any time during stress relaxation after shear cessation, then fracture has 
occurred, otherwise there is no fracture. The velocity profiles shown in Figs. 4 to 6 occur at the time when both /i v+ 
and fi v _ reach their extrema for the case of fracture. When there is no fracture, the velocity profiles are shown when 
either fi v+ reaches its maximum or fi v _ reaches its minimum. When fracture occurs, the position of the fracture plane 
depends on the shape of the specific perturbation. The stress relaxation is independent of the position of the fracture 
plane, as in the experiments of [12] (section III A). 

In about 34% of 300 simulations where A xx , A yy , A xy and 7 are all perturbed simultaneously, the resultant velocity 
profiles resemble the type reported in [12]. The calculations in the manuscript use a set of initial conditions that give 
a fracture with all quantities perturbed, such as subfigure Pi in Fig. 6. 



B. Linear stability analysis 



Linear stability analysis is carried out following the approach described in [21, 29]. To this end the quantities are 
separated into a homogeneous part and a perturbation as 



u(t) =u(t) + 5u(i), 



(11) 



where u(t) = [7, A xx , A xy , A yy ] is a homogeneous state obtained by solving 



d t A xx =2A X!/ 7 - A. 



2r d 



,, - — [1-A] [(/3A + l) A xa: + l] 



2r d 



<9 t A^ =7 + 7 Ayy - A xy - — [1 - A] (/3j4 + 1) A xy 
5 t A ra = - A vv - -± [1 - A\ [{PA + 1) A w + 1] 



A= 1 



TrA 



-1/2 



(12) 



The perturbation 5u(i, y) consists of fluctuations in the velocity gradient direction of the form 



(13) 



Substituting Eq. 8 into Eq. 6 and the momentum equation 



dv 



dt 



+ (vV) 



V • T, 



(14) 



where v is the fluid velocity, p is the fluid density and T is the total stress tensor, gives 
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^A xx (/3 - 1) A 3 - 1 - -^A 3 - 2— [l + (/5 - 1) A] - ^A XX A 4 + 2/3 A 2 
otr Str t r 6 

-k 2 3\ 8A XXik (t) 
+ 2^6A XVtk (t) 

+ b-l)^- 
+ 2A xy S%(t) 



A XX A A 8A XX A 

3t r ' 



-r4 



d t 8A xyik (t) = 



(13-1) 
+ 



-A XV A 8 A A xy 



2/3^-A 

TR 



A 2 - 1 -2^ + 31 -fc 2 P 



(5A„^ fe (t) 



+ 



?+(/?-!) A^A 3 - ^-l3—A xy A i 

otr 6 Tr 



[1 + A yy ]5%(t) 



d t SAyy ik (t) = 



(8 1) ^A yy A 3 ^-A 3 \^AyyX 

OTr 6Tr 6 Tr 



SAyy tk (t) 

SA XXtk (t) 



^{P-l)A y J- 



+28— A 2 - k 2 V 

TR 

k 2 

= - -^8A X y ik (t) - 



A 3 ~ YfAyyA' 

•J'R 3 Tr 

1 - ^-A 3 -2^[l + (8-l)A\- -B^A^A* 

OTR Tr 6 Tr 



<* A TO,fc(*) 



k 2 e 



e = 



~GT d 
9 Cr- 
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where all nonlinear terms have been neglected. In the zero Reynolds number limit p — > this reduces to 



d t 6A X y ik (t) 



d t SAyy tk (t) 



-k 2 V 



1) A - 1 



I^La 3 - 2^ [1 + 09 - 1) A] - l(3A xx A 4 + 2(3 A 2 



SA XXtk (i) 
A 



7' 



e 



(/3 - 1) A ra A 3 - A 3 - ^A ra A 4 



6Ayy tk (t) 



{ p-l)^A xy A 3 - 2 -(3^A 4 A xv 

OT R 3 Tr 



2(3— A 2 - 1 - 2— [1 + (B - 1) A 

TR TR 

?+ 09- 1) t^A^A 3 - ^(3—A xy A 4 
oTr 3 t r 



8A XXtk (t) 
1 + A,, 



5A x?/i fc(t) 



(16) 



SAyy <k (t) 



09-1) 



-A 7 . 7 .A 



-A 



-/3— A yy A 
3 r fl 



<5A rajfe (i) 



^-08- 1) A ra A 3 - 1 - -^A 3 - 2^ [1 + (0 - 1) A] - 2 -(3^-A yy A 4 

OTR 6tr Tr 6 Tr 



-2(3^ A 2 - k 2 V 

TR 



This is a matrix equation of the form 



SAy Vtk (t) 



d t 5u{t) = M(t)-£u(t), 



(17) 



where u = [A^, A xy , A yy \. Similar to the case described in [21], the eigenvalues of the stability matrix M(t) determine 
the (in)stability of the system: the system becomes unstable when the largest real part of an eigenvalue just becomes 
positive [21]. Hence the spinodal (the shear stress at which the fluid goes unstable during startup) for the system can 
be constructed as shown in Fig. 7. This region of instability matches the constitutive curve, similar to the situation 
reported in [21]. 

When the perturbation given in Eq. 9 is used to initialize the system, it induces some inhomogeneity in the system. 
Each point in space can then be considered as a base state and the stability of each of these base states to small 
amplitude fluctuations is also described by the stability matrix M (i) . Hence the most unstable of these base states 
(which is the state whose eigenvalue has the largest real part) can be determined. This approach gives insight into 
the behaviour of the system when the quantities 7, A xx , A xy and A yy are perturbed separately. For 15 different 
initial conditions that give a fracture profile, the eigenvector v m corresponding to the maximum real eigenvalue in 
space at the time of stretch relaxation is heavily dominated by the components A xx and A yy . The components of v m 
for these different initial conditions are shown in Table I, where w^f is the component in the flow direction, is the 
component in the shear direction and v-™ is the component in the velocity gradient direction. Hence perturbing the 
components 7 and A xy separately do not induce 'fracture' (as in Fig. 5) as compared with perturbing the components 
A xx and A yy separately at the same amplitude (as in Fig. 4). 



C. Comparison with experiment 

The calculations in the manuscript are based on the sample SBR 250K whose rheological properties are reported 
in Tables 1 and 2 of [12]. The rheological properties reported in Table 2 of [12] were said to have been measured from 
linear viscoelastic measurements (see section II B of [12]) but the Rouse times reported in Table 2 were estimated 
using Tji — t/(M w /M b ), where r is the reptation time (section II B of [12]). However in our manuscript, the Rouse 
time is calculated using tr = Td/{iZ) (as given in section I of [17]), where Z is the number of entanglements per 
chain. This then implies that the values of rjj quoted in Table 2 of [12] are larger than the values of tr used in our 
manuscript by a factor of 3. We then present the data in [12] into different cases. 

Case I: Intermediate Shear Rate, High Strain- Using (7) = 0.7 s _1 given in Fig. la of [12] and t^ = 4.1s 
quoted in Table 2 of [12] (for the sample SBR 250K) gives (7)7$ ~ 2.9. The sample SBR 250K (see Table 2 of [12]) 
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TABLE I. Components of most unstable eigenvector v m for 15 different initial conditions, for (7} = 200, 70 = 2.5. 



Initial Condition 






«r 


1 


0.9657 


-0.0115 


-0.2594 


2 


0.9767 


-0.0131 


-0.2144 


3 


0.9678 


-0.0118 


-0.2516 


4 


0.9730 


-0.0125 


-0.2304 


5 


0.9632 


-0.0112 


-0.2687 


6 


0.9667 


-0.0117 


-0.2555 


7 


0.9678 


-0.0118 


-0.2513 


8 


0.9636 


-0.0113 


-0.2672 


9 


0.9(39 


-U.Olz ( 


-U.22DD 


10 


0.9683 


-0.0119 


-0.2496 


11 


0.9675 


-0.0118 


-0.2524 


12 


0.9742 


-0.0127 


-0.2255 


13 


0.9666 


-0.0116 


-0.2562 


14 


0.9613 


-0.0110 


-0.2753 


15 


0.9596 


-0.0109 


-0.2812 



has M w = 250000 g/mol and M e = 3300g/mol giving Z = M w /M e = 76. Then using Z = 76, (7)77* = 200 and 
Tci/tr — 3Z gives (7)tr ~ 0.95, which is comparable to the value of (7)77? ~ 1 specihed in case I of the manuscript, 
this is consistent with Fig. 1 of [12]. 

Case II: High Shear Rate, Low Strain- Similarly, (7) = 14 s -1 from [12] gives (7)7^ ~ 57, which is consistent 
with (7)77? > 1 given in case II of the manuscript and it agrees with Fig. 2 of [12]. 

Case III: Low Shear Rate, Low Strain- Again, (7) = 0.05 s _1 gives (7)7-^ ~ 0.2, which is consistent with 
(7)77? < 1 given in case III of the manuscript, this has close agreement with Fig. 7 of [12]. 

The shear stresses at the time of shear cessation for the three cases I, II and III are indicated in Fig. 7. In case I, 
the shear stress had gone through the overshoot and it is beginning to decrease. In case two, the flow is switched off 
before the shear stress reaches the overshoot. Finally in case III, the flow is switched off just before the shear stress 
reaches the overshoot. Figure 1(c) of the manuscript shows a comparison of velocity profiles from the simulations 
and experimental data; the experimental data were obtained from V max in Fig. lc of [12], made dimensionless using 
Vmax = Vmax T /L, where r = 310s (from Table 2 of [12]) and L = 0.7mm as given in section II of [12]. 

Induction time-To check the variation of the delay time after shear cessation before fracture sets in, we performed 
calculations at three different shear rates satisfying (7)77? > 1 (with 77? fixed), similar to Fig. 6 of [12]. For (7) = 600, 
(7)77? ~ 2.8, (7) = 800, (7)773 ~ 3.7 and (7) = 1000, (7)773 ~ 4.6. In all cases, the applied strains are below the strain 
for overshoot at the applied shear rate as indicated by the lines l\ and I2 in Fig. 8(a). The overshoot stress is a linear 
function of the overshoot strain, as in Fig. 6(a) of [12]. Figures 8(bcd) show that, for varying strain and given shear 
rate, the higher plateau stress after stretch relaxation leads to a longer induction time. This characteristic is similar 
to the the situation in the inset of Fig. 6(b) of [12]. 

However, Figs. 8(e-f) show that for fixed strain and varying shear rate, the plateau stresses collapse, and the lower 
applied shear rate yielding a slightly longer induction time for 70 = 2.2. This can be linked to the faster growth rate 
w max observed for the very high shear rates, in which the viscous contribution to the instability dominates. However, 
this behaviour does not match that displayed in the inset of Fig. 6(b) of [12], in which the higher applied shear rate 
resulted in a larger induction time. We do not have an adequate explanation for these discrepancies. 



II. MOVIES 

The movies in https://eudoxus.leeds.ac.uk/dynacop/FracturePage.html illustrate the cases where the fluid 
undergoes fracture after shear cessation (Fracture . avi) and recoil without fracture (Recoil . avi) for case I. To 
achieve this, an initial condition of the form A xx (Q,y) = A(cos(7ry) + 0cos(27ry)) is used to perturb the system. The 
shape and amplitude of this perturbation can be tuned to bring it close to one of the random perturbations which 
yields fracture-like behaviour when the component A xx is perturbed. The amplitude is fixed at A = 0.006 while the 
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FIG. 8. (a) Shear stress versus strain at the three different applied shear rates indicated in the figure such that (7}tr > 1 in all 
cases, red circles: {j)TR — 2.8, green squares: {j)tr — 3.7 and blue diamonds: (7)tr = 4.6. The dashed line connects the strains 
for overshoot and their corresponding stresses for each applied shear rate, while the lines l\ and I2 indicate the applied strains 
70 = 2.2 and 70 = 3.2 respectively, (b)-(f) Stress relaxation after step strains at different applied strains 70 and shear rates (7) 
indicated. Parameters as in Fig. 7. 



parameter <f> is varied to change the shape of the perturbation. The shapes of this perturbation for ip = 0.25 and 
(f> = 0.67 are shown in Fig. 9(a). 

For = 0.67, the fluid fractures after shear cessation, the window on the left of Fracture . avi shows the fluid 
velocity from startup (with the upper plate fixed and the lower plate moving) to shear cessation and continues until 
the end of fracture. Before shear cessation, the fluid is seen to be moving to the left, after which the flow is switched 
off and the velocity vectors go to zero momentarily (except with a slight bulge due to the initial perturbation). The 
sizes of the velocity vectors before shear cessation are larger than their sizes after shear cessation by roughly one order 
of magnitude, hence to make the figure visible in the video, a rescaling of the figure window was carried out after 
shear cessation. The velocity profile v in the video on the left was made dimensionless using v = vTd/L. Then using 
t c i = 310 s and L — 0.7 mm (from [12]) gives the maximum size of velocity vectors v max before shear cessation roughly 
equal to 0.45 mm s -1 and the maximum size after shear cessation is roughly equal to 0.02 mms . The velocity profile 
during fracture is shown in Fig. 9(b). 

The figure window on the right of Fracture.avi shows the corresponding total shear stress T xy /G from startup 
until the end of fracture. The total shear stress builds up quickly when the flow is switched on, and then just after the 
overshoot when the flow is switched off, the total shear stress goes through an initial quick relaxation during which 
the polymer chains relax stretch. It then enters a slow relaxation when reptation sets in. Although some reptation 
had already occurred during stretch relaxation, it becomes the dominant mechanism for stress relaxation after stretch 
relaxation. However, before reptation can completely relax the stress, the growing perturbation causes a sudden quick 
relaxation of stress. By this time the 'fracture plane' is fully developed and the fluid can be seen moving rapidly in 
two different directions on both sides of this plane. Finally when this rapid motion ceases, the stress resumes its slow 
relaxation and the material appears to have healed itself. 

The case of = 0.25, where the peak of the initial perturbation is not sharp enough as in Fig. 9(a), gives a 
completely different relaxation behaviour in the fluid as shown in Recoil.avi. The left window of that figure shows 
the fluid velocity from startup to shear cessation and beyond. Like in the case of (f> = 0.67, the top plate is fixed while 
the lower plate moves to the left. After shear cessation, the perturbation is seen to grow for a while but the fluid 
does not 'break' in two unlike in the case of 4> — 0.67. The growing perturbation loses the competition against the 
background reptation and hence the material heals itself and the fluid velocity vanishes after some time. Like in the 
case of <j> = 0.67, the figure window has been rescaled after shear cessation to make the velocity vectors visible. The 
maximum size of the velocity vectors before shear cessation is roughly equal to 0.45 mms -1 while the maximum size 
after shear cessation is roughly equal to 0.006 mms -1 . The recoil velocity for this case is shown in Fig. 9(b). 
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FIG. 9. (a) Initial perturbation for: <f> = 0.25 and (f> = 0.67, as shown in the movies Recoil.avi and Fracture.avi. (b) Recoil 
after shear cessation for <f> — 0.25 and fracture after shear cessation for <j> — 0.67. Parameters: j3 = 0, Z = Td/(3ri?) = 72, 
(7) = 200, and 70 = 2.5. 



The right window of Recoil.avi shows the corresponding total shear stress for this case. It grows quickly from 
startup like the case of <fr = 0.67, then decays quickly during stretch relaxation and ends up with a slow relaxation due 
to reptation. The stress does not show any stage of rapid relaxation again since reptation is the dominant mechanism 
for stress relaxation in this case. 

The movies were made with a mesh of 100 grid points to reduce the time taken to make them. The relevant 
parameters were set at Z = 72, T> = 10~ 5 , e = 10~ 4 , (3 = 0, 70 = 2.5 and (7) = 200, which represent case I described 
in the manuscript. 
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